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ABSTRACT 

Precise  estimation  of  synchronization  parameters  is  essential  for  re¬ 
liable  data  detection  in  digital  communications  and  phase  errors  can 
result  in  significant  performance  degradation.  The  literature  on  es¬ 
timation  of  synchronization  parameters,  including  the  carrier  fre¬ 
quency  offset,  are  based  on  approximations  or  heuristics  because 
the  optimal  estimation  problem  is  analytically  intractable  for  most 
cases  of  interest.  We  develop  an  online  Bayesian  inference  proce¬ 
dure  for  blind  estimation  of  the  frequency  offset,  for  arbitrary  signal 
constellations.  Our  unified  approach  is  built  on  a  sequential  infer¬ 
ence  procedure  that  leverages  a  novel  result  on  conjugacy  of  the  von 
Mises  and  Gaussian  distributions.  This  conjugacy  allows  for  an  eas¬ 
ily  computable,  closed  form  parametric  expression  for  the  posterior 
distribution  of  the  parameters  given  the  streaming  data,  in  which  hy¬ 
perparameters  are  recursively  updated,  making  the  optimal  sequen¬ 
tial  estimation  problem  mathematically  tractable.  Our  algorithm  is 
computationally  efficient  and  can  be  implemented  in  real-time  with 
very  low  memory  requirements.  Numerical  experiments  are  also 
provided  and  show  that  our  methods  outperform  heuristic  sequen¬ 
tial  carrier  frequency  offset  estimators. 

Index  Terms —  Sequential  Estimation,  Blind  Synchronization, 
Frequency  Offset  Estimation,  Bayesian  Inference,  Phase  tracking. 

1.  INTRODUCTION 

Synchronization  plays  a  critical  role  in  attaining  reliable  digital 
transmission  through  wireless  channels.  Timing  and  frequency  off¬ 
sets  need  to  be  estimated  well  in  order  to  align  the  received  signal 
appropriately  at  the  receiver  before  data  detection.  Optimal  estima¬ 
tors  for  synchronization  parameters  do  not  exist  and  as  a  result,  most 
of  the  existing  literature  uses  approximate  maximum-likelihood 
techniques  and  heuristics  [1,  2],  Even  in  the  simpler  setting,  where 
the  amplitude  is  known  and  constant,  optimal  estimation  of  the 
frequency  offset  is  not  known  in  closed-form  and  only  approxi¬ 
mate  maximum-likelihood  (ML)  and  maximum  a  posterior  (MAP) 
approaches  are  tractable  [3], 

In  digital  communications,  there  is  often  a  mismatch  between 
the  local  oscillator  of  the  transmitter  and  the  receiver.  This  translates 
to  a  carrier  frequency  error  when  downconverting  at  the  receiver, 
which  effectively  rotates  the  signal  constellation  from  sample-to- 
sample.  For  small  carrier  frequency  errors,  A/,  the  constellation 
rotates  slowly  and  tracking  its  rotation  rate  is  required  for  successful 
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data  detection.  Popular  methods  for  carrier  synchronization  include 
the  use  of  a  phase-lock  loop  (PLL).  which  often  requires  hardware 
implementation  or  waveform-level  block-based  processing  in  soft¬ 
ware,  both  of  which  can  be  expensive.  Furthermore,  the  PLL  ap¬ 
proach  often  requires  pilot  data  for  higher  order  modulations.  A 
synchronization  approach  based  on  model-based  sequential  Monte 
Carlo  (SMC)  techniques  was  proposed  in  [4]  to  estimate  the  timing 
offset  and  the  data,  but  is  heavily  dependent  on  known  state-space 
dynamic  models  governing  the  evolution  of  the  parameters,  several 
approximations  are  made  to  make  the  SMC  algorithm  tractable,  and 
is  computationally  expensive  as  the  complexity  of  the  SMC  grows 
exponentially  fast  [4], 

In  contrast,  our  approach  is  not  data-aided  (i.e.,  blind),  works  for 
arbitrary  signal  modulations,  operates  at  the  sample-level,  is  simple 
to  implement,  and  has  very  low  latency  and  storage  requirement  as 
it  proceeds  on-the-fly  in  an  online  fashion.  In  this  paper,  we  pro¬ 
pose  a  novel  sequential  Bayesian  inference  framework  to  estimate 
the  unknown  carrier  frequency  offset  and  the  parameters  of  the  sig¬ 
nal  modulation. 

We  leverage  ideas  from  the  sequential  inference  approach  of  [5] 
to  derive  a  sequential  Bayesian  algorithm  to  estimate  the  frequency 
offset  and  the  constellation  parameters  given  streaming  non-i.i.d.  re¬ 
ceived  samples.  This  is  a  challenging  machine  learning  problem 
primarily  because  of  the  temporal  dynamics  of  the  data  generation 
process.  Our  method  exploits  the  temporal  relation  from  sample- 
to-sample  and  uses  the  online  clustering  and  parameter  estimation 
framework  developed  in  [5].  Solid  empirical  performance  is  ob¬ 
served  for  slow  enough  rotation  rates  in  our  experiments. 


2.  SEQUENTIAL  BAYESIAN  INFERENCE  FRAMEWORK 


Here,  we  review  the  adaptive  sequential  updating  and  greedy  search 
(ASUGS)  framework  of  [5]  for  online  clustering  and  parameter  es¬ 
timation.  Define  the  unknown  parameters  8  =  (/xi , . . . 
where  8  £  [0,  27r)  is  the  unknown  rotation  offset  and  fj,h  is  the 
mean  of  each  class  (i.e..  cluster  center).  Let  the  observations  be 
given  by  y,  £  Rd,  and  7,  to  denote  the  class  label  of  the  *th  ob¬ 
servation  (a  latent  variable).  We  define  the  available  information  at 
time  *  as  yw  =  {yi, . . .  ,y»}  and  7(*_1)  =  {71, . . .  ,7,-1}.  The 
ASUGS  algorithm  is  as  follows.  Set  71  =  1  and  fci  =  1.  Calculate 
7r((?i|yi,7i).  For*  >  2, 

1.  Choose  best  class  label  for  y$: 


7 1 


^2h'  L i ,h' (y i)^ i ,h' 


2.  Update  the  posterior  distribution  using  y i ,  7;: 


7r(07ilyw>7w) «  /(yi|^7i)7r(07ily(l  1))-  (!) 

where  9h  are  the  parameters  of  class  h,  f(yi\0h)  is  the  observa¬ 
tion  density  conditioned  on  class  h.  The  conditional  likelihood 
P(yihi  =  h,  is  denoted  by  Li)h(yi)  and  Tri}h  de¬ 

notes  the  class  priors.  The  algorithm  sequentially  allocates  observa¬ 
tions  y i  to  classes  by  sampling  the  conditional  posterior  probability 
distribution  {q^}- 

3.  PROBLEM  FORMULATION  AND  PROBABILISTIC 
MODEL 


Fig.  1.  Graphical  model  for  simultaneous  phase  tracking  and  param¬ 
eter  estimation  problem. 


Consider  a  digital  communication  system  where  symbols  from  a  dis¬ 
crete  unknown  alphabet,  xm  €  A  are  transmitted  through  a  channel. 
The  received  signal  at  the  front-end  of  the  receiver  consisting  of  a 
matched  filter  can  be  modeled  as: 

y(t)  =  ej2^Aft  ^2  xmg{t  -  mT)  +  w(t) 

m 

where  T  is  the  symbol  period,  g(-)  is  the  raised-cosine  pulse  wave¬ 
form,  and  w(t)  is  additive  white  Gaussian  noise  (AWGN)  with 
power  spectral  density  No/2.  Here,  A /  is  the  carrier  frequency 
error.  Sampling  the  output  of  the  matched  filter  at  a  rate  1  /T,  we 
obtain  the  discrete-time  signal: 

Uk  =  e  xk  +  wk  (2) 

where  yk  =  y(kT),  Wk  =  w(kT),  and  k  =  0, 1, . . .  is  the  discrete¬ 
time  index.  Note  the  phase  rotation  offset  5  =  2tyAJT  is  propor¬ 
tional  to  the  carrier  frequency  error  A/. 

We  consider  a  Bayesian  framework  for  estimating  the  unknown 
mean  for  each  class  and  unknown  phase  rotation  offset  <5.  This 
framework  can  be  extended  to  include  unknown  covariances  as  well 
(corresponding  to  unknown  SNR),  but  is  not  considered  here  for  sim¬ 
plicity  (see  [5]).  The  means  can  be  considered  as  known  by  setting 
the  uncertainty  to  zero  (see  Section  4). 

The  observation  model  (2)  in  vector  form  is  given  by: 


a  normal- von  Mises  joint  distribution.  The  corresponding  graphical 
model  is  shown  in  Fig.  1 . 

For  concreteness,  let  us  write  the  distribution  in  (3): 


f(yi\d)  =p(yi\g,S)  = 


(2na2)dA 


||yi-R(«)Mlli 


p(m) 

p(5) 


1  -^IIm-moII! 

- e  2<to 

(27TCT2)d/2 

1  c^O  cos(<5  — 5q) 

2nI0(K0) 


where  Jo(-)  is  the  modified  Bessel  function  of  order  0. 

Due  to  the  conjugacy  of  the  distributions,  the  posterior  distribu¬ 
tion  7r(0h|y(*_1\  7^1-1^)  always  has  the  form: 

7r(0fc|y(’-1\7(t-1))  =Ar(/rfe|/i^“1),(cr^“1))2I) 

x  V(<S|<S(i-1),K(i-1))  (4) 


where  ^ are  hyperparameters  that  can  be 
recursively  computed  as  new  samples  come  in.  This  factorization  is 
proven  analytically  in  Section  4. 

We  will  show  in  Section  4  that  the  model  (3)  leads  to  closed- 
form  expressions  for  hyperparameter  updates  due  to  conjugacy. 


4.  CONJUGACY  AND  HYPERPARAMETER  UPDATES 


y  i  =  R  (9i)xi  +  w; 

9i+i  =  9i  +  8 


In  this  section,  we  derive  the  hyperparameter  updates.  For  simplicity, 
let  y  be  a  generic  observation  of  the  form: 


where  x,  £  A  are  symbols  from  a  constellation,  w,  is  additive  Gaus¬ 
sian  noise  with  covariance  a2 1.  Here,  R(0)  is  a  rotation  matrix  given 
by 


R(f) 


cos  9  —  sin  9 

sin  9  cos  9 


Without  loss  of  generality,  we  assume  9 1  =  0.  The  unknown  pa¬ 
rameters  in  the  model  are  the  constellation  symbols  x  £  A,  and  the 
phase  offset  8.  Let  K  =  \A\  denote  the  size  of  the  symbol  alphabet. 
The  probabilistic  model  for  the  unknown  parameters  is  given  as: 


y|p,  <5  ~  A/”(jR((5)/r,  cr2I) 

M  ~  A/'(-|/uo,  ctoI) 

<5  ~  V(-|<50,  Ko)  (3) 


where  S)  denote  the  multivariate  normal  distribution  with 

mean  p  and  covariance  matrix  X,  and  V(j<5o,  Ko)  denotes  the  von 
Mises  distribution  with  direction  parameter  <5o  and  concentration  pa¬ 
rameter  Ko-  The  parameters  9  =  ( fj.,8 )  E  Rd  x  [ — 7r,  tt)  follow 


y  =  R(<5)x  +  w 


(5) 


The  conditional  distribution  p(p\8,  y)  oc  p( y\p,  5)p(p)  is  mul¬ 
tivariate  normal  with  mean  and  covariance  given  by: 


E[/r|<5,y]  = 


^%R(5)Ty  +  " 


O-  +  0-,1 


Cov(/x|<5,  y)  =  °  a°  2 1 


Recall  the  von  Mises  prior  distribution  on  8,  i.e.,  p(8 )  oc 
e“°  cosU— ^0)  7^  posterior  is  given  by 


p(8 |y)  oc  p{y\8)p(S) 
oc  exp 


— Ily  -  RWmoIIs 


2  (o'2  +  05) 


+  Kq  COS (J  —  Jo) 


oc  exp 


yTR(<5)/x0 
a2  +  a2 


+  K 0  cos(<5  —  80) 


(6) 


Next,  write  y  =  [ya,  yi]T ,  yo  =  [ Mo, it ,  yo,i]T ■  Then,  direct  calcu¬ 
lations  yield: 

yTH(S)yo  =  ( yTyo )  cos  8  +  (yT yo)  sin  S  (7) 

where  yo  =  [— yo,i,  yo,R]T ■  Furthermore,  from  a  trigonometric 
identity,  we  obtain: 

Ko  cos (8  —  Jo)  =  Ko  cos(Jo)  cos  J  +  k0  sin(Jo)  sin  5  (8) 

Using  (7)  and  (8)  into  (6),  we  obtain: 


p(J|y)  oc  exp  (  (k0  cos(J0)  + 


yTyo 

a2  +  ai 


cos  8 


+  Ko  sin(<50)  + 


T  ~ 

y  mo 
a2  +  aZ 


sin  8 


This  equals  exp(Knew  cos(J  —  J„e,„))  for  all  <5  £  [— n,  n)  iff: 

T 

/  c  \  y  mo  /  c  \ 

K,q  COS  (do)  o  '  o  ^ new  COS(one-U;) 

cr2  +  crg 

T  ~ 

•  /  r-  \  y  Mo  •  /  c  \ 

ttO  Sin(do)  +  ^  j  o"  —  Knew  SlIl(Onet/;) 

CT2  +  C7g 

Solving  for  /tneu,  and  dne™,  we  obtain: 


^nei u  (  ^0  COS (d0)  H- 


T  \  2 

y  mo 


5  new  —  tan 


°  +  °o 
'  k0  sin(J0)  + 

k0  cos  (J0)  + 


+  (  K0  sin(Jo)  +  ^  ~~ 9 
cr2  +  cr^ 


Returning  to  the  model  at  the  ith  time  instant,  pre-multiplying 
with  the  rotation  matrix  R(#i_i)T: 


R(6»i_i)Tyi  =  R(#i-i)TR(0i)xj  +  R(0,_i)tw! 

=  R(6»i_1)-1R(6li_1)R(J)xi  +  R(^_i)tWi 

=  R(J)xi  +  R(#i_i)TWi 


where  R((9i_i)TWi  is  Gaussian  noise  with  zero  mean  and  covari¬ 
ance  <j2I,  due  to  rotation  invariance.  Thus,  the  model  (5)  becomes 
applicable  when  applied  to  R(#,_i)Tyi.  This  essentially  amounts 
to  rotating  the  vector  y,  by  an  angle  9i-\.  A  graphical  illustration 
of  this  rotation  is  shown  in  Fig.  2. 


Fig.  2.  Illustration  of  pre-rotation  of  observation  yi  to  reduce  to 
model  (5). 


To  summarize,  once  the  7ith  component  is  chosen,  the  parame¬ 
ter  updates  for  the  7ith  class  become: 


8 ^  =  tan  1 


<x2+(<4‘1))2 


k(1_1)  cosd(l_1)  +  - - Tx - - 

°‘2  +  (ctt’_  )2  / 


(9) 


(k*-1-1)2  =  ^k^  ^sinj^  1) + 

+  ^K*'1-11  COS  + 

0(i)  =  6|G-1)  +  $(0 


R(6»(l  1))Tyi,/i^  11 

^2  +  (4_1))2 

R(^(i-1))Tyi,4r1) 

*2  +  (4_1))2 


(i)  ( i ) 

r  ~  9),  - 


^hTi,/l(R(0(i))Tyi)7ri,h 


^  +  Ki  ) 


(a«)2  =  (a'r15)2 


<-2  +  (<  J) 


(10) 

(11) 


(i-1) 

(i-i) 


^2  +  (4i 1))2 


(12) 

(13) 


For  exactly  known  means  (i.e.,  constellation  parameters),  one  can 
set  a ^  =  0  and  yj’1  =  yn  for  all  h,  i.  In  this  case,  the  updates 
(12)-(13)  become  superfluous. 


4.1.  Calculation  of  Conditional  Likelihood 

To  calculate  the  class  posteriors  {g^l>},  the  conditional  likelihoods  of 

y  i  =  R(6»(i,)Ty  i  given  assignment  to  class  h  and  the  previous  class 
assignments  need  to  be  calculated  first.  The  conditional  likelihood 
of  y i  given  assignment  to  class  h  and  the  history  (y(i-1l,  7^-1^)  is 
given  by: 

Li,h(yi)  =  j  /(yi|^)7r(6>h|y(l_1),7(l_1))d6»ft  (14) 

We  recall  from  (4)  that  the  posterior  distribution  has  the  product 
form: 

7r(^h|y('_1)i7(i_1))  =  -V (yuly^1^,  (<rjyi))2I)V(J|J(i-1),  k^-1) 

For  large  k^1_1\  which  is  the  case  after  a  few  iterations,  this  can  be 
approximated  as: 

Tr(0h|y(i-1),7(i-1))  *  (<rt”)3)IMS  ~  <5(i_1)) 

(15) 

where  V(-)  denotes  the  Dirac  function.  Using  this  approximation 
(15)  into  (14): 

Li,h{ Yi)  ~  [  A/'(y!!R(0<,,)/i,  cr2I) 

J  Kd 

_  ||R(9W|»yi-^i-1)||| 

=  _ 1 _ g  2(ct2+{,£-1))2) 

(2n(a2  +  (*t1])2))d/2 


(16) 


which  corresponds  to  the  Gaussian  distribution  y i 

(^■1))2)I)- 


5.  SIMULATIONS 

In  this  section,  we  perform  a  simulation  experiment  on  detecting 
and  estimating  the  constellation  parameters  of  a  16-QAM  modula¬ 
tion.  The  phase  rotation  angle  per-sample  is  5  =  3°.  The  correct 
number  of  classes  for  this  data  set  is  16,  each  one  corresponding  to 
4  information  bits. 

Data  symbols  {x;}  are  randomly  (uniformly)  chosen  from  the 
16-QAM  constellation  and  corrupted  with  additive  noise  at  SNR  of 
15dB.  The  data  is  plotted  in  Fig.  3,  along  with  the  original  signaling 
locations  arranged  in  a  rectangular  grid. 
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Fig.  4.  Resulting  compensated  constellation  (left)  and  frequency  off¬ 
set  estimation  performance  (right)  for  naive  sequential  phase  estima¬ 
tor  (17). 
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Fig.  5.  Resulting  compensated  constellation  (left)  and  frequency  off¬ 
set  estimation  performance  (right)  for  sequential  Bayesian  estimator 
(9)-(13). 


Fig.  3.  Noisy  constellation  data  from  16-QAM  constellation  with 
phase  rotation  rate  of  <5  =  3°  at  15dB. 


For  comparison,  we  consider  the  naive  sequential  phase  estima¬ 


tor: 


e(i)  =  z 


Vi  \ 

PA{e-je<-'~1)  yi)  j 


(17) 


where  y j  =  i/f  +jyl,  FU(z)  :=  yh-,  h*  =  argmin^  \\z- yh\\22  is 
the  minimum  norm  projection  on  the  constellation  A.  Fig.  4  shows 
the  clustering  performance  of  this  naive  method,  and  the  clusters  are 
very  noisy  and  not  estimated  well.  As  shown  in  Fig.  4.  the  phase 
offset  S  is  not  learned  as  the  number  of  samples  grow;  although  the 
mean  is  around  3°,  the  variance  is  quite  high  and  is  not  decaying  as 
a  function  of  iteration. 

To  alleviate  this  problem,  we  aim  to  use  the  online  Bayesian  es¬ 
timation  algorithm  developed  in  this  paper  to  estimate  the  frequency 
offset.  Fig.  5  shows  the  clustering  performance  and  the  estimation 
performance  of  the  Bayesian  algorithm.  In  our  simulations,  we  im- 
plemented  the  paremeter  updates  (9)-(l  1),  and  set  y h  =  yh,<rh  = 
0  for  all  h,  i  because  no  uncertainty  in  the  constellation  parameters 
was  assumed.  The  algorithm  was  initialized  with  <5  ■1)  =  0.  With  the 
Bayesian  learning  algorithm  based  on  the  von  Mises  prior,  asymp¬ 
totic  learning  occurs  for  the  phase  offset  S  as  more  samples  are  pro¬ 
cessed.  This  leads  to  a  significantly  more  stable  performance  when 
compared  to  the  naive  phase  estimator  (17). 


6.  CONCLUSION 

We  have  proposed  an  online  Bayesian  framework  for  blind  estima¬ 
tion  of  the  frequency  offset  and  the  parameters  of  an  arbitrary  signal 
constellation.  Our  approach  leverages  novel  conjugate  prior  distri¬ 
bution  theory  for  the  von  Mises  and  Gaussian  distribution,  which 
allows  us  to  derive  closed-form  updates  of  the  hyperparameters  of 
the  posterior  distribution  of  the  unknown  parameters  given  stream¬ 
ing  data.  Simulations  show  that  our  estimation  algorithm  results  in 
fast  convergence  and  learning  of  the  frequency  offset,  and  signifi¬ 
cantly  outperforms  heuristic  online  frequency  offset  estimators. 
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